{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image quilting\n",
    "\n",
    "In this tutorial, we demonstrate the use of the image quilting solver in GeoStats.jl. Image quilting is the fastest multiple-point simulation solver in the literature. For more information about the algorithm and its performance, please watch the following video:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/YJs7jl_Y9yM?rel=0&amp;controls=0&amp;showinfo=0\" frameborder=\"0\" allowfullscreen></iframe>"
      ],
      "text/plain": [
       "HTML{String}(\"<iframe width=\\\"560\\\" height=\\\"315\\\" src=\\\"https://www.youtube.com/embed/YJs7jl_Y9yM?rel=0&amp;controls=0&amp;showinfo=0\\\" frameborder=\\\"0\\\" allowfullscreen></iframe>\")"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using IJulia\n",
    "\n",
    "HTML(\"<iframe width=\\\"560\\\" height=\\\"315\\\" src=\\\"https://www.youtube.com/embed/YJs7jl_Y9yM?rel=0&amp;controls=0&amp;showinfo=0\\\" frameborder=\\\"0\\\" allowfullscreen></iframe>\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we proceed, please install the following packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage GeoStats is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage ImageQuilting is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage GeoStatsImages is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage Plots is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage PyPlot is already installed\n",
      "\u001b[39m"
     ]
    }
   ],
   "source": [
    "for pkg in [\"GeoStats\", \"ImageQuilting\", \"GeoStatsImages\", \"Plots\", \"PyPlot\"]\n",
    "    Pkg.add(pkg)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem definition\n",
    "\n",
    "### Conditional simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2D SimulationProblem (conditional)\n",
       "  data:      4×3 GeoDataFrame (x and y)\n",
       "  domain:    250×250 RegularGrid{Float64,2}\n",
       "  variables: facies (Int64)\n",
       "  N° reals:  3"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using GeoStats\n",
    "\n",
    "# create some artificial data\n",
    "data = DataFrames.DataFrame(x=[50.,190.,150.,150.], y=[50.,50.,170.,190.], facies=[1,0,1,1])\n",
    "\n",
    "geodata = GeoDataFrame(data, [:x,:y])\n",
    "domain  = RegularGrid{Float64}(250, 250)\n",
    "condsim = SimulationProblem(geodata, domain, :facies, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unconditional simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2D SimulationProblem (unconditional)\n",
       "  data:      nothing\n",
       "  domain:    250×250 RegularGrid{Float64,2}\n",
       "  variables: facies (Float64)\n",
       "  N° reals:  3"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uncondsim = SimulationProblem(domain, :facies => Float64, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the problem\n",
    "\n",
    "Like most other multiple-point simulation solvers, image quilting is parametrized with a training image. The [GeoStatsImages.jl](https://github.com/juliohm/GeoStatsImages.jl) package provides various training images from the literature for fast experimentation in Julia. **Please give credit to the data sources if you use these images in your research**.\n",
    "\n",
    "We load a famous image from the geostatistics literature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using GeoStatsImages\n",
    "using Plots; pyplot()\n",
    "\n",
    "TI = training_image(\"Strebelle\")\n",
    "\n",
    "heatmap(TI[:,:,1]', color=:inferno, aspect_ratio=:equal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and define our solver:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using ImageQuilting\n",
    "\n",
    "solver = ImgQuilt(:facies => @NT(TI=TI, template=(30,30,1)));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solver can be used for unconditional simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solution = solve(uncondsim, solver)\n",
    "plot(solution, color=:inferno, size=(1000,300))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "as well as conditional simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solution = solve(condsim, solver)\n",
    "plot(solution, color=:inferno, size=(1000,300))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance issues\n",
    "\n",
    "Currently, there is an implementation detail in ImageQuilting.jl that needs to be fixed in order for the `ImgQuilt` solver in GeoStats.jl to be as performant as the standalone `iqsim` function from that package. For small simulations (e.g. 2D images), the standalone function `iqsim` should perform better and is recommended instead. This issue should be fixed at some point in the future."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
